knitr::opts_chunk$set(eval = FALSE)
Load Libraries
We have three sources, aqautic, terrestrial, and aquatic mosses (a, t, m) and three variables we can measure, C:N, carbon isotopes, and nitrogen isotopes. We want to know, for a given amount of input of a, t, and m, what do we expect our C:N, \(\delta^{13}C\), and \(\delta^{15}N\) to be? Now, since there are errors associated with the proxy values for each input, we’re going to do a Monte Carlo simulation that will select from a distribution within those error ranges and give us a range of outputs.
Define a simple mass balance equation that applies for all proxies.
Mass balance equation:
\[ \begin{align*} x =& \ \textrm{proxy} \ (\delta^{13}C, \ \delta^{15}N, \ or \ C\!:\!N) \\ a =& \ \textrm{amount of algal input} \\ t =& \ \textrm{amount of terretrial plant input} \\ m =& \ \textrm{amount of macrophyte/moss input} \\ x_{total} * (a + t + m) =& (a * x_a + t * x_t + m * x_m) \end{align*} \]
Task 1a: define the normalization function (the evil plankton told me to do it) Task 1b: use the normalization function in the calculate_proxy_error function
#General normalization equation
normalize <- function(x, y, z) {
x / (x + y + z)
} #Write this function!
#General mass balance equation
calculate_proxy <- function(a, proxy_a, t, proxy_t, m, proxy_m) {
(a * proxy_a + t * proxy_t + m * proxy_m) / (a + t + m)
}
#Error for mass balance eqn
calculate_proxy_error <- function(a, proxy_a_error, t, proxy_t_error, m, proxy_m_error) {
a_norm = normalize(a, t, m)
t_norm = normalize(t, a, m)
m_norm = normalize(m, t, a)
sqrt((a_norm * proxy_a_error) ^ 2 + (t_norm * proxy_t_error) ^ 2 + (m_norm * proxy_m_error) ^ 2 )
}
From Florian et. al. (2015)
source_data = read.csv("data/qpt_sources_all.csv")
source_data
Q: Which source should have the biggest (or most unique) affect on C:N? d13C? d15N?
A:
Q: How do the magnitudes of the errors compare across the proxies and accross the sources?
A:
Q: What are the consequences of these uncertainties?
A:
#source parameters
parameters <-
data_frame(
cn_a = 8.8,
cn_t = 48,
cn_m = 31,
d13c_a = -30.7,
d13c_t = -26.9,
d13c_m = -21.8,
d15n_a = 6.5,
d15n_t = -1.6,
d15n_m = 0.77
)
parameter_errors <-
data_frame(
cn_a = 0.4,
cn_t = 29,
cn_m = 4.6,
d13c_a = 0.28,
d13c_t = 1.21,
d13c_m = 1.55,
d15n_a = 0.42,
d15n_t = 3.6,
d15n_m = 0.47
)
Task: Pick values for a, t, and m to model
#Table of input scenarios
input_scenarios <-
expand.grid(
a = seq(0, 100, 10),
t = seq(0, 100, 10),
m = seq(0, 100, 10)
) %>%
merge(parameters) %>%
tbl_df()
#calculate expected proxy values given input scenarios
simulations <-
input_scenarios %>%
merge(data_frame(i = 1:5)) %>%
mutate(
a_norm = normalize(a, t, m),
t_norm = normalize(t, a, m),
m_norm = normalize(m, a, t),
cn_a_w_error = cn_a + rnorm(length(a), sd = 0.25),
cn_t_w_error = cn_t + rnorm(length(t), sd = 15),
cn_m_w_error = cn_m + rnorm(length(m), sd = 4),
cn = calculate_proxy(a, cn_a_w_error, t, cn_t_w_error, m, cn_m_w_error),
d13c_a_w_error = d13c_a + rnorm(length(a), sd = 0.2),
d13c_t_w_error = d13c_t + rnorm(length(t), sd = 0.6),
d13c_m_w_error = d13c_m + rnorm(length(m), sd = 1.5),
d13c = calculate_proxy(a, d13c_a_w_error, t, d13c_t_w_error, m, d13c_m_w_error),
d15n_a_w_error = d15n_a + rnorm(length(a), sd = 0.2),
d15n_t_w_error = d15n_t + rnorm(length(t), sd = 1.8),
d15n_m_w_error = d15n_m + rnorm(length(m), sd = 0.5),
d15n = calculate_proxy(a, d15n_a_w_error, t, d15n_t_w_error, m, d15n_m_w_error)
)
simulations
simulations %>%
ggplot() +
aes(cn, t_norm, color = a_norm) +
geom_point()
Q: Did we simulate the space well enough?
A:
q <- plot_ly(simulations, x = ~cn, y = ~d13c, z = ~d15n, color = ~a_norm) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = "C:N"),
yaxis = list(title = 'delta13C'),
zaxis = list(title = 'delta15N')))
r <- plot_ly(simulations, x = ~cn, y = ~d13c, z = ~d15n, color = ~t_norm) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = "C:N"),
yaxis = list(title = 'delta13C'),
zaxis = list(title = 'delta15N')))
s <- plot_ly(simulations, x = ~cn, y = ~d13c, z = ~d15n, color = ~m_norm) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = "C:N"),
yaxis = list(title = 'delta13C'),
zaxis = list(title = 'delta15N')))
q
r
s
Q: What are we looking at here?
A:
Q: Why does it taper off when it gets yellow?
A:
Now we’re going to work with some actual data that I obtained. The data come from a downcore lake sediment record of lake Qaupat (QPT) on southern Baffin Island, which we looked at on Tuesday. Here’s a map:
knitr::include_graphics("figures/Baffin_Lakes_Map.png")
…and here’s the previous work that’s been done:
Now let’s import the EA data (i.e. C:N, \(\delta^{13}C\) and \(\delta^{15}N\)). I threw in the BSi data from above as well for a comparison.
qpt_data = read.csv("data/QPT16_2A_EA_results_no_outliers.csv")
qpt_data <-
qpt_data %>%
mutate("d15N" = d15N - 25, #This is a relic of the data processing, should fix elsewhere
"Bsi" = parse_number(Bsi)
)
qpt_data
qpt_data_plot <-
qpt_data %>%
gather(var, val, c("Bsi", "d13C", "d15N", "C.N")) %>%
ggplot() +
aes(Age, val, group = var) +
geom_smooth(method = "loess", span = 0.2) +
geom_point() +
facet_grid(var~., scales = "free_y", switch = "y") +
scale_x_reverse(expand = c(0,500)) +
scale_y_continuous(name = "") +
xlab("Age (years)") +
labs(title = "Lake QPT")
qpt_data_plot
Q: Does the EA data align with the BSi? What’s the same and what’s different?
A:
dataplot <-
plot_ly() %>%
add_markers(qpt_data, x = qpt_data$C.N, y = qpt_data$d13C, z = qpt_data$d15N, color = qpt_data$Age) %>%
layout(scene = list(xaxis = list(title = "C:N"),
yaxis = list(title = "delta13C"),
zaxis = list(title = "delta15N")
)) %>%
colorbar(title = "Age")
dataplot
Q: Which variables explain most of the changes? Which explains the least?
A:
Q: Is there a clear age trend? What’s your quick hypothesis as to what’s causing it?
A:
Now let’s pot our data onto our simulation space and see what we find.
q2 <-
plot_ly() %>%
add_markers(simulations, x = simulations$cn, y = simulations$d13c, z = simulations$d15n, color = simulations$a_norm, opacity = 0.3) %>%
add_markers(qpt_data, x = qpt_data$C.N, y = qpt_data$d13C, z = qpt_data$d15N, color = "red") %>%
layout(scene = list(xaxis = list(title = "C:N"),
yaxis = list(title = 'delta13C'),
zaxis = list(title = 'delta15N')))
r2 <-
plot_ly() %>%
add_markers(simulations, x = simulations$cn, y = simulations$d13c, z = simulations$d15n, color = simulations$t_norm, opacity = 0.3) %>%
add_markers(qpt_data, x = qpt_data$C.N, y = qpt_data$d13C, z = qpt_data$d15N, color = "red") %>%
layout(scene = list(xaxis = list(title = "C:N"),
yaxis = list(title = 'delta13C'),
zaxis = list(title = 'delta15N')))
s2 <-
plot_ly() %>%
add_markers(simulations, x = simulations$cn, y = simulations$d13c, z = simulations$d15n, color = simulations$m_norm, opacity = 0.3) %>%
add_markers(qpt_data, x = qpt_data$C.N, y = qpt_data$d13C, z = qpt_data$d15N, color = "red") %>%
layout(scene = list(xaxis = list(title = "C:N"),
yaxis = list(title = 'delta13C'),
zaxis = list(title = 'delta15N')))
q2
r2
s2
Q: Do the data fall within the simulation space?
A:
Q: If no, why might that be the case?
A:
So now we’re going to try to quantify what we just observed visually.
We’re going to perform an optimization test that is going to take a real data point (C:N, d13c, d15n) and find the optimal values of a, t, m (weighted by the errors in the proxies).
How do we know how well our model is doing? We’ll use the formula for distance on a Cartesian space to see how far the optimal points are from the real data. If the distance is zero, the model is right on top of our data point and we have a perfect fit!
\[ distance = \sqrt{(x_{model} - x_{real})^2 + (y_{model} - y_{real})^2 +(z_{model} - z_{real})^2} \]
But being really far off in \(\delta^{13}C\) isn’t so bad because there’s a lot of variation. In other words, being 2 permille off for carbon is not nearly as bad as being 2 permille off for nitrogen.
To deal with this, we’ll also normalize these distances by their standard deviations. It’s a little hand-wavy, but should give us a better answer than no error consideration at all.
\[ \begin{align*} x =& \ \delta^{13}C \\ y =& \ \delta^{15}C \\ z =& \ C:N \\ distance =& \sqrt{(x_{model} - x_{real})^2 / (SD_x)^2 + (y_{model} - y_{real})^2 / (SD_y)^2 +(z_{model} - z_{real})^2 / (SD_z)^2} \end{align*} \]
Here we’ll define our optimization function. It takes a give combination of sources (c(a, t, m)) and a real piece of data (with index = i) and spits out the distance between the value predicted by the sources and the real value.
optim_function <- function(x, i) {
a <- x[1]
t <- x[2]
m <- x[3]
#calculate proxies
dcn <- calculate_proxy(a, parameters$cn_a, t, parameters$cn_t, m, parameters$cn_m)
d13c <- calculate_proxy(a, parameters$d13c_a, t, parameters$d13c_t, m, parameters$d13c_m)
d15n <- calculate_proxy(a, parameters$d15n_a, t, parameters$d15n_t, m, parameters$d15n_m)
#calculate errors
dcn_error <- calculate_proxy_error(a, parameter_errors$cn_a, t, parameter_errors$cn_t, m, parameter_errors$cn_m)
d13c_error <- calculate_proxy_error(a, parameter_errors$d13c_a, t, parameter_errors$d13c_t, m, parameter_errors$d13c_m)
d15n_error <- calculate_proxy_error(a, parameter_errors$d15n_a, t, parameter_errors$d15n_t, m, parameter_errors$d15n_m)
#calculate distance btwn test point and data
sqrt(((dcn - qpt_data$C.N[i]) / dcn_error)^2 + ((d13c - qpt_data$d13C[i]) / d13c_error)^2 + ((d15n - qpt_data$d15N[i]) / d15n_error)^2)
}
Let’s examine what this puts out:
best_fit_i = 3
best_fit <- optim(c(0.33, 0.33, 0.34), i = best_fit_i, optim_function, method = "L-BFGS-B", lower = 0, upper = 1)
best_fit
Now we’re going to take our full data set and find the optimal values of a, t, and m.
qpt_data <- qpt_data %>% mutate(i = row_number())
optim_results <-
qpt_data %>% select(i) %>%
group_by(i) %>%
do({
i <- .$i[1]
result <- optim(c(0.33, 0.33, 0.34), i = i, optim_function, method = "L-BFGS-B", lower = 0, upper = 1)
# normalize output
data_frame(
a_optim = result$par[1] / sum(result$par),
t_optim = result$par[2] / sum(result$par),
m_optim = result$par[3] / sum(result$par),
distance_optim = result$value,
convergence = result$convergence # If != 0, didn't converge
)
})
qpt_data_optim <-
left_join(qpt_data, optim_results, by = "i") %>%
mutate(tot_optim = a_optim + t_optim + m_optim)
qpt_data_optim
Here, we’ll create a ternary diagram with a, t, m (norm) as the axes and all our data points plotted. We’ll color them by age and size them by their best fit distances.
qpt_data_optim_plot <- qpt_data_optim %>%
ggtern(aes(a_optim, t_optim, m_optim, color = Age, size = distance_optim)) +
theme_nomask() +
geom_point() +
scale_color_gradient2(low = "dark blue", mid = "dark green", high = "gold", midpoint = 3500) +
scale_size_continuous(range = c(10,1))
qpt_data_optim_plot
Q: Was our model able to fit the points well? Any points that look suspect?
A:
Q: Do the suspect points make sense? Think back to the 3D plots from earlier.
A:
Now we have some results, but we already suspect that some of the fits we got from our model might not be so great. Let’s analyze how good these fits are.
Task: Let’s see how this works for data points i = 3, 15, and 20 (and 5 if there’s time)
best_fit_i = 20
best_fit <- optim(c(0.33, 0.33, 0.34), i = best_fit_i, optim_function, method = "L-BFGS-B", lower = 0, upper = 1)
#best_fit
best_fit_a = best_fit$par[1] / sum(best_fit$par)
best_fit_t = best_fit$par[2] / sum(best_fit$par)
best_fit_m = best_fit$par[3] / sum(best_fit$par)
This is useful, but it’s still just spitting out a number. It doesn’t really give us a sense of how variable the fit is. So let’s do some plotting.
Now let’s make a ternary plot for (a, t, m), calculate how well each (a, t, m) fits sample i, and display the result as a tile weighted by the goodness of fit (0 = perfect). We’ll also plot our result from “best_fit”, which should show up at the optimal point.
This will allow us to see the gradient of our optim function, i.e. how much wiggle room we have with our fit. It might even allow us to see some second-best fits, i.e. local minima.
Change the resolution of our plot by playing with the number of bins.
expand.grid(
a = seq(0,1, by=0.05),
t = seq(0,1, by=0.05),
m = seq(0,1, by=0.05)
) %>%
tbl_df() %>%
mutate(
allsum = a+t+m,
a = a/allsum,
t = t/allsum,
m = m/allsum,
value = map2_dbl(a, t, ~optim_function(c(.x, .y, best_fit_m), best_fit_i))) %>%
filter(!is.na(value), !is.nan(value)) %>%
ggtern(aes(a,t,m)) +
# calculate mean estimate in each bin
geom_tri_tern(bins = 5, mapping = aes(fill=..stat..,value=value), fun = mean) +
# mark minimum
geom_point(data = qpt_data_optim, aes(
x = qpt_data_optim$a_optim[best_fit_i],
y = qpt_data_optim$t_optim[best_fit_i],
z = qpt_data_optim$m_optim[best_fit_i]), color = "red")
# + scale_fill_gradientn(colours = rainbow(7))
Q: What did you learn about the goodness of our fits?
A:
Q: Did it look like all of the points had a global minimum on our simulation space?
A:
Q: What conclusions can you draw? How well did this model work?
A:
Outside of this file, I also played a lot with a mixing model called MixSIAR. We’ll discuss it briefly.
MixSIAR is a Bayesian Markov Chain Monte Carlo stable isotope mixing model. Like the simulation we just did, it’s going to take us from our data to our a, t, m space, but in a more mathematically refined way.
First, we’ll make a plot analogous to the one in the man-eating lions paper (Yeakel et. al., 2009):
isospace <-
ggplot() +
geom_point(data = qpt_data, aes(x = d13C, y = d15N, color = Age)) +
#geom_path(data = qpt_data, aes(x = d13C, y = d15N, color = Age), size = 0.2) +
scale_color_gradient2(low = "dark blue", mid = "dark green", high = "gold", midpoint = 3500) +
geom_point(data = source_data, aes(x = Meand13C, y = Meand15N), color = "black") +
geom_errorbar(data = source_data, aes(
x = Meand13C, y = Meand15N, ymin = Meand15N - SDd15N, ymax = Meand15N + SDd15N)) +
geom_errorbarh(data = source_data, aes(
x = Meand13C, y = Meand15N, xmin = Meand13C - SDd13C, xmax = Meand13C + SDd13C)) +
geom_text(data = source_data, aes(x = Meand13C, y = Meand15N, label = Source, hjust = -0.5, vjust = 2)) +
labs(
title = "Lake QPT Isospace",
color = "Age (yrs)",
x = latex2exp::TeX("$\\delta^{13}C$"),
y = latex2exp::TeX("$\\delta^{15}N$"))
isospace
Q: Spoiler alert: the model didn’t really work out. Why might that be?
A:
Previously, it looked like C:N changed more than \(\delta^{15}N\). Let’s see if that still looks true on an isoscape:
#source_data_cn <- read.csv("data/qpt_sources_cn.csv")
isospace_cn <-
ggplot() +
geom_point(data = qpt_data, aes(x = d13C, y = C.N, color = Age)) +
#geom_path(data = qpt_data, aes(x = d13C, y = d15N, color = Age), size = 0.2) +
scale_color_gradient2(low = "dark blue", mid = "dark green", high = "gold", midpoint = 3500) +
geom_point(data = source_data, aes(x = Meand13C, y = MeanCN), color = "black") +
geom_errorbar(data = source_data, aes(
x = Meand13C, y = MeanCN, ymin = MeanCN - SDCN, ymax = MeanCN + SDCN)) +
geom_errorbarh(data = source_data, aes(
x = Meand13C, y = MeanCN, xmin = Meand13C - SDd13C, xmax = Meand13C + SDd13C)) +
geom_text(data = source_data, aes(x = Meand13C, y = MeanCN, label = Source, hjust = -0.5, vjust = 2)) +
labs(
title = "Lake QPT Isospace",
color = "Age (yrs)",
x = latex2exp::TeX("$\\delta^{13}C$"),
y = latex2exp::TeX("$\\C:N$"))
isospace_cn
Q: Does this look better? Think we’d get a good model output from this?
A: